MODULE md1_hermite_mpoles

   USE md1_global_consts
   USE md1_global_types
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: md1_build_1D_mpoles

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE md1_build_1D_mpoles(xyz,batch,origin,Qjie)
   
      USE md1_hermite_expansion, ONLY: md1_build_Ecoefs

      IMPLICIT NONE
      INTEGER(INTK),        INTENT(IN)  :: xyz
      TYPE(md1_prim_batch), INTENT(IN)  :: batch(:)
      REAL(REALK),          INTENT(IN)  :: origin(3)
      REAL(REALK),          INTENT(OUT) :: Qjie(:,0:,0:,0:)
      
      REAL(REALK), ALLOCATABLE :: Mte(:,:,:)
      REAL(REALK), ALLOCATABLE :: Etij(:,:,:,:)
      INTEGER(INTK) :: np, Iang, Jang, LMAX

      np   = UBOUND(Qjie,1)
      Jang = UBOUND(Qjie,2)
      Iang = UBOUND(Qjie,3)
      LMAX = UBOUND(Qjie,4)

      ALLOCATE (Mte(np, -1:LMAX+1, -1:LMAX+1))
      ALLOCATE (Etij(np, 0:Iang+Jang, 0:Iang, 0:Jang))

      CALL md1_build_hermite_mpoles(batch,origin,xyz,LMAX,Mte)
!if (xyz == 2) then
!print *, "MteY"
!do i = 1, np
!   print *, i
!   print *, Mte(i,0:,0:)
!   print *
!end do
!end if
      CALL md1_build_Ecoefs(batch,xyz,Iang,Jang,Etij)
      CALL md1_transform_hermite_to_car(LMAX,Etij,Mte,Qjie(:,:,:,:))
!if (xyz == 2) then
!print *, "QjieY"
!do i = 1, np
!   print *, i
!   print *, Qjie(i,0:,0:,0:)
!   print *
!end do
!end if

      DEALLOCATE (Mte)
      DEALLOCATE (Etij)

   END SUBROUTINE md1_build_1D_mpoles

!-------------------------------------------------------------------------------

   SUBROUTINE md1_build_hermite_mpoles(batch,origin,xyz,LMAX,Mte)

      IMPLICIT NONE
      TYPE(md1_prim_batch), INTENT(IN)  :: batch(:)
      REAL(REALK),          INTENT(IN)  :: origin(3)
      INTEGER(INTK),        INTENT(IN)  :: xyz, LMAX
      REAL(REALK),          INTENT(OUT) :: Mte(:,-1:,-1:)

      INTEGER(INTK) :: e, t
      REAL(REALK) :: expI(SIZE(batch))
      REAL(REALK) :: PC(SIZE(batch))

      ! Set multipole origin
      PC(:) = batch(:)%P(xyz) - origin(xyz)

      expI(:) = half / batch(:)%exp
      Mte(:,:,:) = zero
      Mte(:,0,0) = SQRT(two*PI*expI(:))

      DO e = 1, LMAX
         DO t = 0, e
            Mte(:,t,e) =  PC(:) * Mte(:,t,e-1)           &
                      + expI(:) * Mte(:,t+1,e-1)         &
                      + DBLE(t) * Mte(:,t-1,e-1)
         END DO
      END DO

   END SUBROUTINE md1_build_hermite_mpoles

!-------------------------------------------------------------------------------

   SUBROUTINE md1_transform_hermite_to_car(LMAX,Etij,Mte,Qjie)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: LMAX
      REAL(REALK),   INTENT(IN)  :: Etij(:,0:,0:,0:)
      REAL(REALK),   INTENT(IN)  :: Mte(:,-1:,-1:)
      REAL(REALK),   INTENT(OUT) :: Qjie(:,0:,0:,0:)

      INTEGER(INTK) :: i,j,e,t,tmax

      Qjie(:,:,:,:) = zero

      DO e = 0, LMAX
         DO i = 0, UBOUND(Etij,3)
            DO j = 0, UBOUND(Etij,4)
               tmax = MIN((i+j),e)
               DO t = 0, tmax
                  Qjie(:,j,i,e) = Qjie(:,j,i,e) + Etij(:,t,i,j) * Mte(:,t,e)
               END DO
            END DO
         END DO
      END DO

   END SUBROUTINE md1_transform_hermite_to_car

!-------------------------------------------------------------------------------

END MODULE md1_hermite_mpoles

